function prof_data = profjoin(prof_1, prof_2, draw_index, num_pts);

% PROFJOIN - joins separate profiles
%
%    PROF = PROFJOIN(PROF1, PROF2, INDEX, PTS) joins two sets of profiles
%    (PROF1 & PROF2) to form the output (PROF).  INDEX contains the indices
%    of the profile drawings in the global structure DRAW_STRUCT and PTS
%    determines the number of points used to in the overlap to anneal the
%    profiles smoothly.

% By:   S.C. Molitor (smolitor@med.unc.edu)
% Date: August 18, 2000

% initialize output

prof_data = [];

% validate arguments

if (nargin ~= 4)
   msgbox('Invalid number of arguments', 'PROFJOIN Error', 'warn');
   return
elseif (~isnumeric(prof_1) | isempty(prof_1))
   msgbox('PROF1 must be a numeric array', 'PROFJOIN Error', 'warn');
   return
elseif (~isnumeric(prof_2) | isempty(prof_2))
   msgbox('PROF2 must be a numeric array', 'PROFJOIN Error', 'warn');
   return
elseif (size(prof_1, 2) ~= size(prof_2, 2))
   msgbox('PROF1 and PROF2 must have the same number of columns', 'PROFJOIN Error', 'warn');
   return
elseif (~isnumeric(draw_index) | (length(draw_index) ~= 2))
   msgbox('INDEX must be a two-element vector', 'PROFJOIN Error', 'warn');
   return
elseif (~isscalar(num_pts))
   msgbox('PTS must be a scalar', 'warn');
   return
end

% obtain info from DRAW_STRUCT
% make sure specified drawings are profiles

global DRAW_STRUCT
if (length(DRAW_STRUCT) < max(draw_index))
   msgbox(sprintf('DRAW_STRUCT has %d drawings', length(DRAW_STRUCT)), 'PROFJOIN Error', 'warn');
   return
elseif (~isfield(DRAW_STRUCT, 'type') | ~isfield(DRAW_STRUCT, 'shape'))
   msgbox('DRAW_STRUCT missing TYPE/SHAPE fields', 'PROFJOIN Error', 'warn');
   return
elseif (~isfield(DRAW_STRUCT, 'x') | ~isfield(DRAW_STRUCT, 'y'))
   msgbox('DRAW_STRUCT missing X/Y fields', 'PROFJOIN Error', 'warn');
   return
elseif (~strcmp(DRAW_STRUCT(draw_index(1)).type, 'PROFILE'))
   msgbox(sprintf('Drawing %d is %s', draw_index(1), DRAW_STRUCT(draw_index(1)).type), 'PROFJOIN Error', 'warn');
   return
elseif (~strcmp(DRAW_STRUCT(draw_index(2)).type, 'PROFILE'))
   msgbox(sprintf('Drawing %d is %s', draw_index(2), DRAW_STRUCT(draw_index(2)).type), 'PROFJOIN Error', 'warn');
   return
end

% obtain coordinates from each drawing
% find minimum distance between coordinates

xy_1 = double(drawmask(DRAW_STRUCT(draw_index(1)).type, ...
   DRAW_STRUCT(draw_index(1)).shape, DRAW_STRUCT(draw_index(1)).x, DRAW_STRUCT(draw_index(1)).y));
xy_2 = double(drawmask(DRAW_STRUCT(draw_index(2)).type, ...
   DRAW_STRUCT(draw_index(2)).shape, DRAW_STRUCT(draw_index(2)).x, DRAW_STRUCT(draw_index(2)).y));
[x1, x2] = meshgrid(xy_1(1, :), xy_2(1, :));
[y1, y2] = meshgrid(xy_1(2, :), xy_2(2, :));
dist2 = (x2 - x1).^2 + (y2 - y1).^2;
min_index = find(dist2 == min(min(dist2)));
[min_index(2), min_index(1)] = ind2sub(size(dist2), min_index(1));

% create arrays to indicate indices for overlap & separation
% prevent array indices from being out of bounds

overlap(1, :) = [ceil(-num_pts/2) floor(num_pts/2)] + min_index(1);
overlap(2, :) = [ceil(-num_pts/2) floor(num_pts/2)] + min_index(2);
if ((min(overlap(1, :)) <= 1) | (max(overlap(1, :)) > size(prof_1, 1)))
   msgbox('Overlap indices for PROF1 out of bounds.  Try reducing PTS.', 'PROFJOIN Error', 'warn');
   return
elseif ((min(overlap(2, :)) < 1) | (max(overlap(2, :)) >= size(prof_2, 1)))
   msgbox('Overlap indices for PROF2 out of bounds.  Try reducing PTS.', 'PROFJOIN Error', 'warn');
   return
elseif (diff(overlap(1, :)) ~= diff(overlap(2, :)))
   msgbox('PROF1 & PROF2 must have the same number of overlap indices', 'PROFJOIN Error', 'warn');
   return
end
alone(1, :) = [1 (min(overlap(1, :)) - 1)];
alone(2, :) = [(max(overlap(2, :)) + 1) size(prof_2, 1)];

% convert first columns to distance increment
% will be cumulatively summed after profiles are joined

prof_1(:, 1) = [0; diff(prof_1(:, 1))];
prof_2(:, 1) = [0; diff(prof_2(:, 1))];

% pad profile arrays with zeroes to align

add_1 = [prof_1(alone(1, 1) : overlap(1, 2), :); zeros(diff(alone(2, :)) + 1, size(prof_1, 2))];
add_2 = [zeros(diff(alone(1, :)) + 1, size(prof_2, 2)); prof_2(overlap(2, 1) : alone(2, 2), :)];

% create weighting array
% add profiles as weighted sum
% convert first column back to cumulative distance

weight = [zeros(diff(alone(1, :)) + 1, 1); ...
      ([0 : diff(overlap(1, :))]' + 1) / (diff(overlap(1, :)) + 1); ones(diff(alone(2, :)) + 1, 1)];
for i = 1 : size(prof_1, 2)
   if (i == 1)
      prof_data = cumsum((1 - weight) .* add_1(:, i) + weight .* add_2(:, i));
   else
      prof_data(:, i) = (1 - weight) .* add_1(:, i) + weight .* add_2(:, i);
   end
end
return
